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If the fermion mass is large enough, the phase of the fermion determinant of QCD at finite density is strongly 
correlated with the imaginary part of the Polyakov loop. This fact can be exploited to reduce the fluctuations 
of the phase significantly , making numerical simulations feasible in regions of parameters where the naive brute 
force method does not work. 



1. The Polyakov loop and the phase of the 
determinant 

It is well known that the euclidean path in- 
tegral representation of the QCD partition func- 
tion at finite chemical potential suffers from the 
so-called sign problem: the fermion determinant 
is complex and the theory cannot be simulated 
with the usual Monte Carlo method. One could 
still try to get results with the simple brute 
force method, that is, simulate a positive mea- 
sure which contains the pure gauge action and 
the modulus of the determinant, treating the 
phase as an observable. Then, the expectation 
value of any observable O is given by the ratio 
(O) = (O expi8) m / (cos 8) m , where (•),„ denotes 
the expectation value using the modulus of the 
determinant as a probability measure, and 9 is the 
phase of the determinant (PD). Unfortunately, 
the expectation value of cos 6 is a positive quan- 
tity exponentially small with the volume. Since 
the relative error (O) is given by 



8(Q) _ S{Oexpi6) m 5(cos9) r 
W <0expi0) m + (cos0}„. 
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(cos#) m must be measured very accurately to 
achieve a given accuracy for (0). This requires 
statistic growing exponentially with the volume. 
Obviously, this is, from the numerical point of 
view, an almost hopeless task. 

Due to the up to now unsurmountable difficul- 
ties with QCD with light quarks at finite density, 
increasing attention is being paid to the limit of 
infinitely heavy quarks The sign problem 



still remains in this limit, but numerical compu- 
tations are easier For large quark mass m, 
the logarithm of the fermion determinant can be 
expanded in powers of l/(2m) (we focus our dis- 
cussion on staggered fermions). The first term 
sensitive to the chemical potential \x is 



e 
2m 



where Nt is the number of lattice points in the 
temporal direction, x labels the sites in a given 
temporal slice, and is the ordered product 
of all temporal links attached to x. Hence, to 
this order, the PD is proportional to the imag- 
inary part of the Polyakov loop (IPPL): 9 = 
cV s Li, with Li — l/V s ^2 s ImTiLg and c = 
2 sinh(^iYr)/(2TO) ArT . Of course, higher order 
corrections will destroy this linear relationship. 

The static limit (m, /x — ■> oo with c = 
[exp((i)/2m] NT fixed) has been studied in An- 
alyzing the data of pj , we found a strong correla- 
tion between the PD and the IPPL, as can be seen 
in Figs. 1 and 2 0). The strong correlation dis- 
played in Fig. 1 is not surprising, since the data 
correspond to c w 0.05, which is small enough for 
the linear relation previously discussed to become 
essentially exact. Fig. 2 is more interesting, since 
c ~ 1.66 is not small. There, we can see that the 
linear correlation still holds, though the width of 
the band is much larger than in Fig. 1. Very re- 
cently, a paper confirming these findings from a 
continuum analysis has appeared p] . 
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Figure 1. Correlation between The PD 9 and The 
IPPL Li in static QCD, with C = 1/20 and /3 = 
5.0 on a lattice L 3 x 2. The data are from Ref. 1. 

2. Fixing the imaginary part of the 
Polyakov loop 

Given the correlation between the PD and the 
IPPL, it is plausible that, at least for heavy 
enough quarks, the fluctuations of the PD would 
be suppressed to a large extent provided we con- 
strain our path integral to configurations with 
real Polyakov loop. To see whether this is pos- 
sible we write the partition function as Z = 
J dpi exp (—V s J 7 (pi)), where V s is the spatial lat- 
tice volume and 

e~ v ^ = J [dU]e- s *M det A(U) 5{ Pl - U) . (3) 

For V s — ► oo the integral in pi is saturated by 
the saddle point p° p , which is in general complex, 
since T is. 

It is easy to show that the saddle point is indeed 
purely imaginary. The expectation value of the 
IPPL coincides with p\ p , and therefore 

sp _ {Li cos 6) m . {Li sm0) m 
Pl ~ {cos9) m +t {cos0) m • ( ' 

For each gauge configuration we also have its com- 
plex conjugate, for which Li changes sign. From 
the loop expansion of the logarithm of the fermion 
determinant, we see that the modulus and the 
phase depend only on the real and imaginary part 
of the loops respectively. Hence, the modulus 
and the pure gauge action do not change, while 
8 changes sign. The first term in the r.h.s. of (Q) 
vanishes, so that pl F is purely imaginary. 
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Figure 2. Same as Fig. 1 with C « 1.66 on a 
6 3 x 4 lattice. 



The existence of a saddle point for the partition 
function implies the equivalence between canoni- 
cal and microcanonical ensembles. One can con- 
strain an observable to its saddle point value. 
This only makes sense when the saddle point is 
real. Therefore, we can only constrain the IPPL 
to zero in those cases where its expectation value 
is zero. We expect this to happen at zero tem- 
perature [0. At finite temperature, the expec- 
tation value of the Polyakov loop gives the free 
energy of a heavy quark, and its complex conju- 
gate that of a heavy antiquark. Since these two 
free energies should be different, the expectation 
value of the IPPL cannot be zero. However, if 
the temperature is small, the expectation value 
of IPPL should be exponentially small with 1/T. 
The static limit at strong coupling, which can be 
solved analytically jjj, confirms this. In this case, 
pl" = ( c 2 - c)/(c 3 + 1), with c = {e»/2m) NT , in 
agreement with our expectation. Therefore, we 
can constrain the IPPL to zero in simulations of 
the low temperature phase of QCD at finite den- 
sity. 

3. Constrained Monte Carlo 

To test our ideas we developed a HMC algo- 
rithm in which the molecular dynamics is forced 
to evolve on the surface of zero IPPL. This is 
easily achieved by introducing a Lagrange multi- 
plier which must be computed at each step of the 
molecular dynamics by solving a non-linear equa- 



3 



tion, which is the condition that the constraint 
must be obeyed at each step. It can be shown 
that the dynamics is reversible and that detailed 
balance is satisfied. The additional cost in CPU 
time caused by the constraint is small: 10% more 
than the unconstrained case for quenched simu- 
lations, and completely negligible with dynamical 
fcrmions. 

As a preliminary test, we made several 
quenched runs to check the gain of a constrained 
simulation in comparison with an unconstrained 
one. We worked with a 4 3 x 6 lattice at = 1.0 to 
ensure that we were in the low temperature phase. 
We got a set of 300 equilibrium quenched config- 
urations and we diagonalized its associated mass- 
less staggered fermion matrix, for \x — 0.5, 1.0, 1.5 
and 2.0. Having the eigenvalues of the massless 
matrix, it is very cheap to get the determinant 
for any value of the mass. 

Let us describe our results. If the mass is large, 
there is no sign problem in the small and large [i 
regions. Simulations are easy but not interesting 
there. When /j, is of the order m the sign problem 
becomes severe: this is the region where the onset 
transition, separating the zero from the finite den- 
sity phases, takes place. It is very interesting to 
determine it accurately. Fig. 3a displays (cos 9) 
as a function of m, for /i = 1.5. The logarithmic 
scale allows us to see the relative error entering 
Eq. (|l|). In some cases, the relative error in the 
unconstrained simulation is one order of magni- 
tude larger than in the constrained one. Notice 
the severe sign problem signaling a transition for 
m ~ 1.75. The transition region is very broad in 
the unconstrained simulation, and in fact covers 
most of the finite density region, between satu- 
ration (low m) and zero density (high m). Ac- 
tually, almost nothing can be inferred about the 
finite density phase in the unconstrained case. In 
the constrained simulation, the sign problem oc- 
curs in a much narrower window, so that the on- 
set transition as well as the finite density phase 
could be analysed. From Fig. 3a, we see that 
fi = 1.5 is nearly critical for m S [1.75, 1.82]. As 
expected, the constraint in the IPPL improves 
the sign problem and makes numerical simula- 
tions feasible. Fig. 3b displays the distribution 
of 9/tt for jjL — 1.5 and m — 1.5. The difference 
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Figure 3. a) Average of cos 9 vs. fermion mass 
for n = 1.5; b) Histogram of 9/tt for [i = 1.5 and 
m = 1.5 . 



between the constrained and unconstrained cases 
is manifest. An extended version of this work has 
recently appeared 

We thank the authors of Ref. Q|, especially 
Doug Toussaint, for making their data available 
to us. Ph. de F. thanks Mike Creutz for helpful 
discussions. V.L. acknowledges useful discussions 
with R. Aloisio, V. Azcoiti and A. Galante. 

REFERENCES 



T.C. Blum, J.E. Hetrick, and D. Toussaint, 

Phys. Rev. Lett. 76 (1996) 1019. 

J. Engels et al., hep-lat/9903030 and these 

proceedings. 

R. Aloisio et al. Phys. Lett. B453 (1999) 275; 
hep-lat/9903004 . 

Ph. de Forcrand and V. Laliena, [hep- 



lat/9907004 



5. K. Langfeld, G. Shin, hep-lat/9907006 



